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ABSTRACT 

A  least-squares  technique  for  the  estimation  of 
the  phase  velocity  and  apparent  azimuth  of  a  plane  wave 
propagating  across  a  horizontal  array  is  developed  and 
applied  to  two  events:  the  Rayleigh  waves  from  an 
Alaskan  earthquake  recorded  on  the  seven-element  IJ BO 
long-period  seismic  array,  and  the  acoustic-gravity  waves 
from  a  presumed  atmospheric  nuclear  explosion, 
recorded  on  part  of  the  thirteen-element  microharo- 
graphic  array  at  LAMA0  Both  time-domain  and  frequency- 
domain  formulations  of  the  technique  were  used,  i0e„, 
we  used  both  the  multichannel  correlation  and  spec¬ 
tral  matrices  to  measure  the  interchannel  time  shifts  . 

The  results  suggest  that  the  time-domain  application 
may  provide  a  useful  algorithm  for  signal  detection. 

The  phase  velocity  dispersion  results  in  both  cases 
agree  reasonably  well  with  theory  and  with  group 
velocities  calculated  by  another  method.  Prewhitening 
(flattening  the  spectral  peaks  by  linear  filtering) 
was  tested  as  a  means  of  preventing  spectral  leakage 
from  contaminating  the  derived  phase  velocity  and 
azimuth;  the  results  suggest  that  prewhitening  decreases 
the  standard  error  when  the  s ignal - to-noise  ratio  is 


high,  but  is  ineffective  when  the  signal-to-noise  rati 
is  low. 
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1.1  ST  01’  FIGURES 


Figure  Title  Figure  No. 

Notation  and  sign  conventions  for  a 

plane  wave  propagating  across  an  array.  1 

UBO  recording  of  Rayleigh  waves  from  the 

Alaskan  earthquake  of  2  June  1969  (see 

Table  1)  .  2 

LAMA  recording  of  acoustic-gravity  waves 
from  a  presumed  atomspheric 

explosion  (see  Table  1).  - 

Signal - to-noise  ratio  for  the  Alaskan  earth¬ 
quake  Rayleigh  waves,  with  an  estimate  of 
group  velocity  added.  4 

Arrival  azimuth  dispersion  for  the  Alaskan 
earthquake  Rayleigh  waves  (no  prewnitening) .  5 

Signal- to-noisje  ratio  for  the  Alaskan  earth¬ 
quake  Rayleigh  waves.  6 

Phase  velocity  dispersion  for  the 
acoustic-gravity  waves  (no 

prewhitening) .  7 

Arrival  azimuth  dispersion  for  the 
acoustic-gravity  waves  (no 

prewhitening) .  8 

Time-varying  spectrum  ("seismoprint")  of 

the  Alaskan  earthquake  Rayleigh  waves.  9 

Time-varying  spectrum  of  the 

acoustic-gravity  waves.  10 
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INTRODUCTION 

McCowan  and  Flinn  (1968)  described  a  least-squares 
method  for  estimating  the  N-l  interchannel  time  shifts 
between  corresponding  phases  of  a  plane  wavefront  cross¬ 
ing  a  horizontal  array  of  N  elements.  The  technique 
made  use  of  estimates  of  the  interchannel  time  shifts 
based  on  the  time  lags  between  all  possible  pairs  of 
elements  as  measured  from  the  multichannel  cross¬ 
correlation  matrix.  It  was  shown  that  this  led  to  N/2 
degrees  of  freedom  for  estimation  of  each  of  the  time 
shifts . 

In  this  paper  we  discuss  u  simple  extension  of  this 
method:  the  direct  least-squares  estimation  of  the 
propagation  velocity  and  apparent  arrival  azimuth  of 
the  hypothetical  propagating  wavefront.  Since  in  this 
case  there  are  only  two  variables  to  solve  for  (instead 

of  N-l),  the  number  of  degrees  of  freedom  is  increased 

2 

to  (N  -N-8)/4,  and  thus  we  would  expect  greater  sta¬ 
bility  in  the  estimates,  provided  of  course  that  what 
we  are  observing  is  indeed  a  plane  wave) . 

Another  motivation  for  this  extension  of  the  method 
is  detection  analysis.  At  least  one  detection  algorithm 
now  in  use  (the  N-4  multichannel  correlator;  see  Brown 
1962)  directly  measures  the  signal  phase  velocity  and 
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azimuth  of  arrival.  The  detection  process  includes  measure¬ 
ment  of  the  persistency  in  time  of  the  propagation  para¬ 
meters.  The  least-squares  method  presented  here  has 
similar  characteristics,  but  in  addition  to  the  estima¬ 
tion  of  the  propagation  parameters  themselves,  the 
present  technique  gives  estimates  of  the  standard  errors 
of  the  derived  variables.  We  might  therefore  expect 
that  these  standard  errors,  displayed  as  a  function  of 
time,  may  also  be  a  practical  detection  algorithm,  since 
presumably  they  should  \e  significantly  smaller  when  a 
signal  is  present  than  when  the  records  contain  random 
and  incoherent  noise. 

In  this  paper  we  describe  and  demonstrate  this 
method  of  estimating  the  propagation  parameters  and 
their  errors.  A  more  detailed  discussion  and  additional 
results  are  given  by  McCowan  and  Flinn  (1970)., 
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THEORY 

When  a  plane  wavefront  crosses  a  horizontal  array 
of  N  elements,  as  illustrated  in  Figure  1,  the  time 
interval  between  the  arrival  of  a  corresponding  wave 
feature  at  two  different  elements  is: 


«ij *  <yv  ■  <*i-y  y-y  cose 


v 


(l) 


Since  this  equation  must  hold  between  all  possible  pairs 
of  elements,  the  problem  of  calculating  V  and  9  from  the 
observed  is  clearly  over-determined,,  For  example,  with 
an  array  of  four  elements,  the  equations  are: 


where  the  unknowns  are: 

f  -  S  COS0 

1  '  V  ■»  f2  =  ~~v 
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PROPAGATION 

DIRECTION 


ELEMENT  i 


AYCOS0 


ELEMENT  j 


Figure  1.  Notation  and  sign  conventions  for  a  plane  wave  propagat 
across  an  array. 


2 

There  are  (N  -N)/2  equations  when  there  are  N  array 

2 

elements,  so  there  are  (N  -N-8)/4  decrees  of  freedom 
left  over  to  provide  information  about  the  errors  in 
the  estimates  z^  and  f2  (Anderson,  1958). 

In  matrix  form,  equation  (2)  is: 

Hi  -  t 

and  the  least-squares  estimate  of  the  unknown  vector  is 


(4) 


~ 

Providing  the  errors  are  small,  reasonable  estimates  of 
the  propagation  phase  velocity  and  azimuth  can  be  obtained 
from  the  relations: 


(5) 


/N  S\  r\  Ss  f\  y 

v  =  Cf t  +  fir* 


-i  /N  /N 

0  =  tan_1(f1/f2) 


(6) 


Again  from  least-squares  theory,  the  variances  and 

A  A 

covariance  in  the  original  unknowns  f^  and  f2  are: 


2  1 


o/=  i  [t't  -  f'  (H'H)f][(H'K)"3  ]  _ 

±  U  **  *•  ~  I  J  1 


(7) 


o22=  i  [t»t-ff (H»H)f )[(H»H)  132  2 


(8) 
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Thus  equations  (6)  and  (12)  -  (15)  constitute  a  least- 

squares  method  Tor  obtaining  estimates  of  the  phase 

velocity  and  arrival  azimuth  of  a  propagating  waveform, 

together  with  their  respective  variances,  provided  we 

2 

can  reliably  measure  the  (N  /N)/z  interelement  time 

delays,,  If  the  errors  are  small,  the  linear  Taylor  series 

expansion  is  valid  out  to  several  standard  deviations „ 

2 

Thus  if  t  has  a  normal  distribution,  so  do  V  and  0, 
approximately  o 

There  arc  two  direct  ways  to  measure  these  delays. 
The  firsf  is  to  compute  the  cross-correlation  matrix 
and  measure  the  delays  by  finding  the  time  lag  where 
each  cross -cor relation  function  is  a  maximum.  Thi' 
method  has  the  advantage  that  the  averaging  implicit 
in  computing  the  correlation  functions  smoothes  out 
the  effect  of  uncorrelated  noise.  Measuring  the  delays 
in  this  way  can  be  expected  to  give  values  of  phase 
velocity  and  arrival  azimuth  which  are  averages  ov-  r 
all  frequencies. 

The  other  way  to  measure  the  delays  is  to  compute 
them  from  the  cross -spectral  matrix 

,i, 

S.  .  (w)  =  E[X*.'(u)X.  ( U) )  3 
1]  1  3 


— —  1 1  til  Ml 


wheie  I:  [  ]  denotes  the  expected  value  and  the  asterisk 
indicates  complex  conjugation.  The  spectrum  of  a  propa¬ 
gating  signal  is: 

sij(‘*))  =  aO)  exp  [  iwx  A  j  (oi)  ]  (IS) 

where  x^  is  the  time  delay  between  array  elements  i  and  j . 

One  dilficulxy  is  the  necessity  to  make  sure  that  the 
phase  spectrum  for  each  cross-spectrum  is  properly 
unwound",  i.e.,  the  periodicity  of  phase  removed.  This 
can  lead  to  trouble  when  two  elements  of  the  array  are 
far  apart.  If  difficulties  arise  for  this  reason,  however, 
then  che  time-domain  calculation  is  also  ambiguous, 
because  the  cross-correlation  functions  will  be  longer 
and  there  will  be  an  ambiguity  of  peak  correlation. 

I  he  major  advantage  of  the  frequency-domain  approach  is 
that  independent  information  is  obtained  at  each  fre¬ 
quency,  provided  of  course  that  the  cross-spectral 
estimates  are  themselves  uncorrelated  between  frequencies. 

A 

To  summarize,  if  the  errors  in  t  are  small  and 
normally  distributed,  we  expect  that  the  estimates  of 
phase  velocity  and  azimuth  of  the  propagating  wavefront 
to  be  approximately  normally  distributed  variables. 
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Therefore  the  standard  deviations  given  in  equations  (12) 
and  (13)  should  represent  familiar  confidence  intervals. 

The  velocity  dispersion  that  is  measured  by  the 
least-squares  technique  is  an  array -d ependent  phenomenon. 
Neither  spectral  matrices  nor  correlation  matrices  con¬ 
tain  any  information  about  path-dependent  dispersion: 
they  are  both  sensitive  only  to  phase  di f fcrcnccs ,  not 
absolute  phase.  Therefore  the  dispersion  results  obtained 
by  the  technique  described  here  depend  only  on  structure 
of  the  propagation  medium  in  the  vicinity  of  the  array. 
This  is  advantageous  in  that  a  localized  measurement  of 
structure  can  be  made,  but  the  advantage  may  be  offset 
by  the  problem  of  resolution:  small  arrays,  though 
yielding  consistent  results,  see  only  a  small  portion 
of  the  wave  propagation  process.  Thus  we  must  weigh  the 
results  of  the  technique  against  the  limited  scope  of 
observation  inherent  in  small  arrays. 

Smart  (1971)  has  shown  that  leakage  of  power  between 
frequencies  can  cause  misinterpretation  of  frequency- 
wavenumber  spectra.  That  is,  an  apparent  energy  peak  on 
the  usual  display  of  power  contoured  as  a  function  of 
two  wavenumber  components  may  not  represent  a  parcel 
of  energy  travelling  with  phase  velocity  equal  to  the 
wavenumber  location  of  the  peak  divided  by  the  frequency 


of  the  cross-section;  the  peak  may  simply  he  leakage  along 
lines  of  constant  wavenumber  from  a  true  spectral  peak 
at  some  other  frequency,,  Smart  lias  proposed  a  method  lor 
eliminating  this  effect0 

We  were  not  sure  that  this  leakage  might  not  cause 
errors  in  the  values  of  phase  velocity  and  azimuth  derived 
by  the  method  presented  in  this  report,  since  the  spectra 
of  the  Rayleigh  and  acoustic-gravity  waves  do  contain 
large  peaks.  To  test  the  leakage  offer  we  performed 
all  analyses  twice,  once  witli  the  origi  a'l  records,  and 
once  after  "prewhitening",  i„e„,  linear  filtering  to 
lower  the  spectral  peaks  (McGowan  and  Flinn  (1970)  for 
details) . 
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rw°  events  were  selected  to  demonstrate  the  method: 
the  Rayleigh  waves  from  an  Alaska  earthquake  recorded  in 
!,tah  b>'  the  H'lLta  Basin  Observatory  array  of  long-period 
vertical  seismometers,  and  an  acoustic-gravity  wave 
•signal  from  a  presumed  atmospheric  nuclear  explosion,  recorded 
in  Montana  by  part  of  the  barge  Aperture  Microbarograph 
Array  (Cook  and  Bedard,  1971)  These  two  events  were 
chosen  because  of  easy  availability  of  digital  data, 
and  good  signal-to-noise  ratio. 

Information  about  these  events  was  taken  from  the 
NOaA/NOS  Preliminary  Determination  of  Hpicenters  releases; 
a  summary  is  given  in  Table  I.  Relative  locations  of  the 
array  elements  can  be  obtained  from  Carpenter  (1965), 

I  he  records  used  are  shown  in  Figures  2  and  3,  along 
with  the  array  beams  directed  toward  the  NOS  epicenters „ 

LAMA  elements  D1  -  D4  were  not  used  in  this  study., 

The  responses  of  the  two  arrays  are  roughly  circular, 
and  the  diameter  of  the  1  dB  contour  of  both  is  about 
0o09  cycles  per  kilometer. 

The  data  were  detrended  and  the  microbarograph 
records  were  bandpass - f i 1 te red  and  decimated  by  a 
lactor  of  four.  For  the  leakage  effect  experiment, 
the  data  channels  were  prewhitened,  each  independently. 
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For  both  events  the  prewhitening  filter  was  a  fifty-point 
least-squares  approximate  inverse  filtei ,  applied  in  the 
time  domain  before  computing  correlations  or  spectra 
(fcr  background,  see  for  example  Treitel  and  Robinson, 
1966) o  The  cross-correlation  matrices  were  computed  using 
an  algorithm  described  by  Flinn  (1968),  and  the  cross- 
spectral  matrices  using  a  Fourier  transform  method 
described  by  McCowan  (1968) . 

We  estimated  group  velocity  by  two  methods.  First, 
we  used  a  process  described  by  Archambeau  and  Flinn 
(1965)  in  which  the  records  are  put  through  a  comb  of 
very  narrow  recursive  bandpass  filters  and  the  group 
velocity  calculated  from  the  arrival  time  of  maximum 
envelope  amplitude  on  the  filter  outputs.  The  program 
compensates  for  instrument  group  delay,  although  for 
the  epicentral  distances  involved  here  the  effect  of 
instrument  group  delay  is  negligible. 

The  second  method  used  to  determine  group  velocity 
was  to  compute  power  spectra  of  overlapped  time  windows 
and  to  display  the  results  as  a  function  of  frequency 
and  time,  i.e.,  frequency  and  group  velocity.  To  the 
best  of  our  knowledge,  this  technique  to  determine 
dispersion  of  acoustic-gravity  or  seismic  waves 
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RESULTS 


The  results  of  this  study  are  presented  in  Table  JI 
and  Figures  4  through  10 .  For  both  events  the  time-domain 
calculation  (Table  II)  shows  clearly  the  effects  of 
prewhitening,  which  gives  worse  results  for  the  acoustic- 
gravity  waves  and  better  results  for  the  Rayleigh  waves, 
ln  terms  of  standard  errors  of  the  estimates.  To  see 
why  this  is  so  we  must  look  at  the  results  in  the 
frequency  domain. 

The  frequency-domain  results  for  both  events  show 
little  difference  between  the  unprewhitened  and  pre¬ 
whitened  cases  (the  prewhitened  results  are  not  given 
here;  they  are  presented  in  detail  by  Flinn  and  McCowan, 
1970),  The  complexity  and  variability  of  the  phase 
velocity  and  arrival  azimuth  for  the  Rayleigh  waves 
from  the  Alaskan  earthquake  (Figures  4  and  5)  are  remark¬ 
able,  It  is  obvious  that  scatter  of  at  least  one  standard 
deviation  is  present  in  these  data,  since  no  reasonable 
pi ane- layered  earth  structure  yields  a  phase  velocity 
dispersion  curve  given  by  the  mean  of  the  estimates  in 
the  anomalous  region  around  0.06  -  0,07  cps  (which  shows 
up  in  both  the  phase  velocity  and  azimuth  plots).  The 
s ignal- to-noise  ratio  (Figure  6)  appears  to  be  reason¬ 
ably  high,  so  noise  contamination  is  probably  not  the 
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cause  of  the  variation  in  phase  velocity  (sec  below  for 
a  discussion  of  the  group  velocity  data  shown  in  Figure 
4) o  It  is  likely  that  multipath  Rayleigh  wave  arrivals 
aie  present.  These  might  be  caused,  by  lateral  refrac¬ 
tions,  the  effect  of  which  for  dispersed  waves  is  fre¬ 
quency  dependent „  The  difference  in  arrival  azimuth 
between  8  and  30  seconds  period  is  almost  30°.  This  can 
be  explained  by  the  complicated  travel  path  from  Alaska 
to  UBO,  which  roughly  goes  down  the  west  side  of  the 
Rocky  Mountains  and  crosses  them  at  a  highly  oblique 
angle,  which  would  tend  to  magnify  any  refraction 
effects.  Another  obvious  possible  cause  is  derivation 
from  horizontal  layering  along  the  propagation  path„ 

In  contrast,  che  acoustic-gravity  waves  arc 
comparatively  smooth  and  well-behaved  (Figure  7),  and 
the  calculated  arrival  azimuth  is  practically  constant 
over  a  broad  frequency  range  (Figure  8).  It  seems 
therefore  that  these  records  are  relatively  uncon¬ 
taminated  by  winds,  which  would  cause  a  frequency- 
dependent  refraction  similar  to  that  for  the  Rayleigh 
waves.  The  flatness  is  probably  one  reason  why  pre¬ 
whitening  harms  these  results;  another  is  the  com¬ 
paratively  low  signal-to-noise  ratio  of  microbarograph 
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records).  We  notice  from  Table  II  that  the  arrival 
azimuths  for  both  the  prewhitened  and  unprewhitened  cases 
are  quite  close  to  the  great-circle  azimuth.  We  conclude 
that  the  addition  of  white  noise  to  an  event  of  low 
s ignal- to-noise  ratio  should  detract  from  the  accuracy 
of  the  estimates,  and  probably  cannot  enhance  them.  For 
the  Rayleigh  waves,  however,  prewhitening  tends  to  smooth 
out  in  frequency  the  already  high  signal-to-noise  ratio, 
thereby  giving  more  uniform  results. 

In  the  time  domain  the  Rayleigh  wave  phase  velocity 
derived  from  the  unprewhitened  data  is  clearly  too  high, 
but  that  derived  from  the  prewhitened  data  is  not  only 
reasonable  but  agrees  with  the  average  over  frequency  of 
the  frequency-domain  calculation. 

Figures  (4)  and  (7)  also  show  group  velocity  dis¬ 
persion  curves  calculated  by  program  SWAP  (Archambeau  and 
Flinn,  1965).  It  is  interesting  to  note  that  the  pro¬ 
nounced  group  velocity  minima  predicted  by  Press  and 
Harkrider  (1962)  are  clearly  evident  in  the  dispersion 
curve  for  the  acoustic-gravity  waves. 

On  the  basis  of  the  average  slope  of  the  phase 
velocity  curve  we  have  also  calculated  the  group 
velocity  dispersion  for  the  Rayleigh  waves  from  the 
Alaskan  earthquake  (Figure  4).  Not  only  is  the  general 


shape  reasonable,  but  the  actual  values  agree  quite  well 
with  the  group  velocity  calculated  by  SWAP.  In  these 
comparisions  it  must  be  remembered  that  program  SWAP 
measures  dispersion  cause'1  by  group  delay  at  the  source, 
along  the  propagation  path,  and  at  the  receiver,  while 
the  least-squares  technique  presented  here  measures  only 
the  last  effect. 

Figures  9  and  10  show  time-varying  spectra  of  the 
data.  Group  velocities  derived  by  tracing  ridge  lines 
(Figure  4)  agree  reasonably  well  with  those  calculated 
by  SWAP  except  in  the  vicinity  of  the  80-second  and  the 
150-second  periods  for  the  acoustic-gravity  waves,  where 
the  SWAP  group  velocity  is  about  two  percent  higher  than 
the  seismoprint  ridge  line  indicates.  This  is  probably 
due  to  mode  interference,  i.e.,  to  energy  interference 
in  the  narrow-band  filter  outputs,  to  which  the  SWA]’ 

procedure  is  known  to  be  sensitive  (see  Archambeau  and 
Flinn,  1965). 

Other  interesting  features  appearing  on  the  acoustic- 
gravity  wave  seismoprints  are  discussed  by  Smart  and 
Flinn  (1971). 
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CONCLUSIONS 


1°  Frequency-domain  application  of  the  least-squares 
estimation  of  signal  propagation  parameters  can  be  a 
useful  method  for  measuring  phase  velocity  dispersion  from 
array  data.  In  the  case  of  highly  complicated  dispersion 
curves,  such  as  those  observed  for  the  Rayleigh  waves 
studied  here,  multipath  arrivals  may  be  a  significant 
cause  of  scatter. 

2.  The  time-domain  application  is  more  accurate 
when  the  signal - to-noise  ratio  is  high  and  the  disper¬ 
sion  curves  are  relatively  smooth.  This  procedure,  when 
supplemented  by  recursively  computed  correlation  functions 
may  be  the  more  useful  in  signal  detection. 

3.  The  frequency-domain  dispersion  curves  seem  to 
be  affected  by  the  limits  of  the  array  response  at  the 
low-frequency  end  and  the  instrument  system  response  at 
the  high-frequency  end.  This  might  serve  as  a  measure 
of  array  resolution  in  addition  to  the  more  standard 


methods 
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